﻿using System;
using System.Diagnostics.Contracts;

namespace MathService
{
    public static class SpecialFunctions
    {
        //private const double e = Consts.e;
        private const double pi = Consts.pi;
        //private const double pi2 = Consts.pi2;
        /// <summary>Pi/2</summary>
        private const double pi_2 = pi / 2;

        /// <summary>Точность алгоритмов 5e-16</summary>
        public const double Eps = 5E-16;

        private const double c_MaxRealNumber = 1E300;
        private const double c_MinRealNumber = 1E-300;
        private const double c_LogMinRealNumber = 690.77552789821368;
        private const double c_LogMaxRealNumber = -690.77552789821368;
        private const double sqrt2 = 1.4142135623730951;

        /// <summary>Гамма-функция</summary>
        public static class Gamma
        {

            [Pure, System.Diagnostics.DebuggerStepThrough]
            public static double G(double x)
            {
                double z;
                var sgngam = 1;
                var q = Math.Abs(x);
                if(q > 33)
                {
                    if(x >= 0) z = GStir(x);
                    else
                    {
                        double p = (int)Math.Floor(q);
                        var i = (int)Math.Round(p);
                        if(i % 2 == 0) sgngam = -1;
                        z = q - p;
                        if(z > .5) z = q - ++p;
                        z = q * Math.Sin(pi * z);
                        z = Math.Abs(z);
                        z = Math.PI / (z * GStir(q));
                    }
                    return sgngam * z;
                }
                z = 1;

                while(x >= 3) z *= --x;

                while(x < 0)
                    if(x <= -.000000001) z /= x++;
                    else
                        return z / ((1 + .5772156649015329 * x) * x);

                while(x < 2)
                    if(x >= .000000001) z /= x++;
                    else
                        return z / ((1 + .5772156649015329 * x) * x);

                if(Math.Abs(x - 2) < Eps) return z;
                x -= 2;

                var pp = 1.60119522476751861407E-4;
                pp = 1.19135147006586384913E-3 + x * pp;
                pp = 1.04213797561761569935E-2 + x * pp;
                pp = 4.76367800457137231464E-2 + x * pp;
                pp = 2.07448227648435975150E-1 + x * pp;
                pp = 4.94214826801497100753E-1 + x * pp;
                pp = 9.99999999999999996796E-1 + x * pp;

                var qq = -2.31581873324120129819E-5;
                qq = 5.39605580493303397842E-4 + x * qq;
                qq = -4.45641913851797240494E-3 + x * qq;
                qq = 1.18139785222060435552E-2 + x * qq;
                qq = 3.58236398605498653373E-2 + x * qq;
                qq = -2.34591795718243348568E-1 + x * qq;
                qq = 7.14304917030273074085E-2 + x * qq;
                qq = 1.00000000000000000320 + x * qq;

                return z * pp / qq;
            }

            [System.Diagnostics.DebuggerStepThrough]
            private static double GStir(double x)
            {
                var w = 1 / x;
                var stir = 7.87311395793093628397E-4;
                stir = -2.29549961613378126380E-4 + w * stir;
                stir = -2.68132617805781232825E-3 + w * stir;
                stir = 3.47222221605458667310E-3 + w * stir;
                stir = 8.33333333333482257126E-2 + w * stir;
                w = 1 + w * stir;
                var y = Math.Exp(x);

                if(x <= 143.01608) y = Math.Pow(x, x - .5) / y;
                else
                {
                    var v = Math.Pow(x, .5 * x - .25);
                    y = v * v / y;
                }
                return 2.50662827463100050242 * y * w;
            }

            [System.Diagnostics.DebuggerStepThrough]
            public static double LnG(double x)
            {
                int sign;
                return LnG(x, out sign);
            }

            [System.Diagnostics.DebuggerStepThrough]
            public static double LnG(double x, out int sign)
            {
                double p;
                double q;
                double z;
                sign = 1;
                const double logpi = 1.14472988584940017414;
                const double lc_Ls2Pi = .91893853320467274178;
                if(x < -34)
                {
                    q = -x;
                    int tmp_sign;
                    var w = LnG(q, out tmp_sign);
                    p = (int)Math.Floor(q);
                    var i = (int)Math.Round(p);
                    sign = i % 2 == 0 ? -1 : 1;
                    z = q - p;
                    if(z > .5) z = ++p - q;
                    return logpi - Math.Log(q * Math.Sin(pi * z)) - w;
                }
                if(x < 13)
                {
                    z = 1;
                    p = 0;
                    var u = x;

                    while(u >= 3)
                    {
                        u = x + --p;
                        z *= u;
                    }

                    while(u < 2)
                    {
                        z /= u;
                        u = x + ++p;
                    }

                    if(z >= 0) sign = 1;
                    else
                    {
                        sign = -1;
                        z = -z;
                    }

                    if(Math.Abs(u - 2) < Eps) return Math.Log(z);

                    p -= 2;
                    x += p;

                    var b = -1378.25152569120859100;
                    b = -38801.6315134637840924 + x * b;
                    b = -331612.992738871184744 + x * b;
                    b = -1162370.97492762307383 + x * b;
                    b = -1721737.00820839662146 + x * b;
                    b = -853555.664245765465627 + x * b;

                    var c = 1.0;
                    c = -351.815701436523470549 + x * c;
                    c = -17064.2106651881159223 + x * c;
                    c = -220528.590553854454839 + x * c;
                    c = -1139334.44367982507207 + x * c;
                    c = -2532523.07177582951285 + x * c;
                    c = -2018891.41433532773231 + x * c;

                    return Math.Log(z) + x * b / c;
                }
                q = (x - .5) * Math.Log(x) - x + lc_Ls2Pi;
                if(x > 100000000) return q;

                p = 1 / (x * x);
                if(x >= 1000)
                    q += ((7.9365079365079365079365 * .0001 * p - 2.7777777777777777777778 * .001) * p
                        + .0833333333333333333333) / x;
                else
                {
                    var a = 8.11614167470508450300 * .0001;
                    a = -(5.95061904284301438324 * .0001) + p * a;
                    a = 7.93650340457716943945 * .0001 + p * a;
                    a = -(2.77777777730099687205 * .001) + p * a;
                    a = 8.33333333333331927722 * .01 + p * a;
                    q += a / x;
                }
                return q;
            }
        }

        /// <summary>Класс функций Бесселя</summary>
        [Copyright("1984, 1987, 1989, 2000 by Stephen L. Moshier")]
        public static class BesselFunction
        {
            private static void BesselAsympt0(double x, out double P, out double Q)
            {
                var xsq = 64 / (x * x);
                var p2 = .0;
                p2 = 2485.271928957404011288128951 + xsq * p2;
                p2 = 153982.6532623911470917825993 + xsq * p2;
                p2 = 2016135.283049983642487182349 + xsq * p2;
                p2 = 8413041.456550439208464315611 + xsq * p2;
                p2 = 12332384.76817638145232406055 + xsq * p2;
                p2 = 5393485.083869438325262122897 + xsq * p2;
                var q2 = 1.0;
                q2 = 2615.700736920839685159081813 + xsq * q2;
                q2 = 156001.7276940030940592769933 + xsq * q2;
                q2 = 2025066.801570134013891035236 + xsq * q2;
                q2 = 8426449.050629797331554404810 + xsq * q2;
                q2 = 12338310.22786324960844856182 + xsq * q2;
                q2 = 5393485.083869438325560444960 + xsq * q2;
                var p3 = -.0;
                p3 = -4.887199395841261531199129300 + xsq * p3;
                p3 = -226.2630641933704113967255053 + xsq * p3;
                p3 = -2365.956170779108192723612816 + xsq * p3;
                p3 = -8239.066313485606568803548860 + xsq * p3;
                p3 = -10381.41698748464093880530341 + xsq * p3;
                p3 = -3984.617357595222463506790588 + xsq * p3;
                var q3 = 1.0;
                q3 = 408.7714673983499223402830260 + xsq * q3;
                q3 = 15704.89191515395519392882766 + xsq * q3;
                q3 = 156021.3206679291652539287109 + xsq * q3;
                q3 = 533291.3634216897168722255057 + xsq * q3;
                q3 = 666745.4239319826986004038103 + xsq * q3;
                q3 = 255015.5108860942382983170882 + xsq * q3;
                P = p2 / q2;
                Q = 8 * p3 / q3 / x;
            }

            private static void BesselAsympt1(double x, out double P, out double Q)
            {
                var xsq = 64 / (x * x);
                var p2 = -1611.616644324610116477412898;
                p2 = -109824.0554345934672737413139 + xsq * p2;
                p2 = -1523529.351181137383255105722 + xsq * p2;
                p2 = -6603373.248364939109255245434 + xsq * p2;
                p2 = -9942246.505077641195658377899 + xsq * p2;
                p2 = -4435757.816794127857114720794 + xsq * p2;

                var q2 = 1.0;
                q2 = -1455.009440190496182453565068 + xsq * q2;
                q2 = -107263.8599110382011903063867 + xsq * q2;
                q2 = -1511809.506634160881644546358 + xsq * q2;
                q2 = -6585339.479723087072826915069 + xsq * q2;
                q2 = -9934124.389934585658967556309 + xsq * q2;
                q2 = -4435757.816794127856828016962 + xsq * q2;

                var p3 = 35.26513384663603218592175580;
                p3 = 1706.375429020768002061283546 + xsq * p3;
                p3 = 18494.26287322386679652009819 + xsq * p3;
                p3 = 66178.83658127083517939992166 + xsq * p3;
                p3 = 85145.16067533570196555001171 + xsq * p3;
                p3 = 33220.91340985722351859704442 + xsq * p3;

                var q3 = 1.0;
                q3 = 863.8367769604990967475517183 + xsq * q3;
                q3 = 37890.22974577220264142952256 + xsq * q3;
                q3 = 400294.4358226697511708610813 + xsq * q3;
                q3 = 1419460.669603720892855755253 + xsq * q3;
                q3 = 1819458.042243997298924553839 + xsq * q3;
                q3 = 708712.8194102874357377502472 + xsq * q3;

                P = p2 / q2;
                Q = 8 * p3 / q3 / x;
            }

            private static void BesselMFirstCheb(double c, ref double b0, ref double b1, ref double b2) { b0 = c; b1 = 0; b2 = 0; }

            private static void BesselMnExtCheb(double x, double c, ref double b0, ref double b1, ref double b2)
            {
                b2 = b1;
                b1 = b0;
                b0 = x * b1 - b2 + c;
            }

            private static void BesselM1FirstCheb(double c, ref double b0, ref double b1, ref double b2) { b0 = c; b1 = 0; b2 = 0; }

            private static void BesselM1NextCheb(double x, double c, ref double b0, ref double b1, ref double b2)
            {
                b2 = b1;
                b1 = b0;
                b0 = x * b1 - b2 + c;
            }

            /// <summary>Функция Бесселя 0 порядка</summary>
            /// <param name="x">Аргумент</param>
            /// <returns>Значение функции Бесселя нулевого порядка</returns>
            [System.Diagnostics.DebuggerStepThrough]
            public static double J0(double x)
            {
                x = Math.Abs(x);

                if(x > 8)
                {
                    double p0, q0;
                    BesselAsympt0(x, out p0, out q0);
                    var nn = x - pi / 4;
                    return Math.Sqrt(2 / pi / x) * (p0 * Math.Cos(nn) - q0 * Math.Sin(nn));
                }

                var xsq = x * x;

                var p1 = 26857.86856980014981415848441;
                p1 = -40504123.71833132706360663322 + xsq * p1;
                p1 = 25071582855.36881945555156435 + xsq * p1;
                p1 = -8085222034853.793871199468171 + xsq * p1;
                p1 = 1434354939140344.111664316553 + xsq * p1;
                p1 = -136762035308817138.6865416609 + xsq * p1;
                p1 = 6382059341072356562.289432465 + xsq * p1;
                p1 = -117915762910761053603.8440800 + xsq * p1;
                p1 = 493378725179413356181.6813446 + xsq * p1;

                var q1 = 1.0;
                q1 = 1363.063652328970604442810507 + xsq * q1;
                q1 = 1114636.098462985378182402543 + xsq * q1;
                q1 = 669998767.2982239671814028660 + xsq * q1;
                q1 = 312304311494.1213172572469442 + xsq * q1;
                q1 = 112775673967979.8507056031594 + xsq * q1;
                q1 = 30246356167094626.98627330784 + xsq * q1;
                q1 = 5428918384092285160.200195092 + xsq * q1;
                q1 = 493378725179413356211.3278438 + xsq * q1;

                return p1 / q1;
            }

            /// <summary>Функция Бесселя 1 порядка</summary>
            /// <param name="x">Аргумент функции</param>
            /// <returns>Значение функции Бесселя первого порядка</returns>
            [System.Diagnostics.DebuggerStepThrough]
            public static double J1(double x)
            {
                double s = Math.Sign(x);

                x = Math.Abs(x);

                if(x > 8)
                {
                    double p0, q0;
                    BesselAsympt1(x, out p0, out q0);
                    var nn = x - 3 * Math.PI / 4;
                    var result = Math.Sqrt(2 / pi / x) * (p0 * Math.Cos(nn) - q0 * Math.Sin(nn));
                    if(s < 0) result = -result;
                    return result;
                }

                var xsq = x * x;

                var p1 = 2701.122710892323414856790990;
                p1 = -4695753.530642995859767162166 + xsq * p1;
                p1 = 3413234182.301700539091292655 + xsq * p1;
                p1 = -1322983480332.126453125473247 + xsq * p1;
                p1 = 290879526383477.5409737601689 + xsq * p1;
                p1 = -35888175699101060.50743641413 + xsq * p1;
                p1 = 2316433580634002297.931815435 + xsq * p1;
                p1 = -66721065689249162980.20941484 + xsq * p1;
                p1 = 581199354001606143928.050809 + xsq * p1;

                var q1 = 1.0;
                q1 = 1606.931573481487801970916749 + xsq * q1;
                q1 = 1501793.594998585505921097578 + xsq * q1;
                q1 = 1013863514.358673989967045588 + xsq * q1;
                q1 = 524371026216.7649715406728642 + xsq * q1;
                q1 = 208166122130760.7351240184229 + xsq * q1;
                q1 = 60920613989175217.46105196863 + xsq * q1;
                q1 = 11857707121903209998.37113348 + xsq * q1;
                q1 = 1162398708003212287858.529400 + xsq * q1;

                return s * x * p1 / q1;
            }

            /// <summary>Функция Бесселя n порядка</summary>
            /// <param name="n">Порядок функции Бесселя</param>
            /// <param name="x">Аргумент Функции Бесселя</param>
            /// <returns>Значение функции Бесселя n порядка</returns>
            [System.Diagnostics.DebuggerStepThrough]
            public static double Jn(int n, double x)
            {
                int sg;

                if(n >= 0) sg = 1;
                else
                {
                    n = -n;
                    sg = n % 2 == 0 ? 1 : -1;
                }

                if(x < 0)
                {
                    if(n % 2 != 0) sg = -sg;
                    x = -x;
                }

                if(n == 0) return sg * J0(x);

                if(n == 1) return sg * J1(x);

                if(n == 2) return Math.Abs(x - 0) < Eps ? 0 : sg * 2 * J1(x) / x - J0(x);

                if(x < Eps) return 0;

                var k = 53;
                var pk = 2d * (n + k);
                var ans = pk;
                var xk = x * x;

                do
                {
                    pk -= 2;
                    ans = pk - xk / ans;
                    k--;
                } while(k != 0);

                ans = x / ans;
                pk = 1;
                var pkm1 = 1 / ans;
                k = n - 1;
                var r = 2d * k;

                do
                {
                    var pkm2 = (pkm1 * r - pk * x) / x;
                    pk = pkm1;
                    pkm1 = pkm2;
                    r -= 2;
                    k--;
                } while(k != 0);

                return sg * (Math.Abs(pk) > Math.Abs(pkm1) ? J1(x) / pk : J0(x) / pkm1);
            }

            /// <summary>Функция Бесселя второго типа, нулевого порядка</summary>
            /// <param name="x">Аргумент функции Бесселя второго типа, нулевого порядка</param>
            /// <returns>Значение функции Бесселя второго типа, нулевого порядка</returns>
            [System.Diagnostics.DebuggerStepThrough]
            public static double Y0(double x)
            {
                if(x > 8)
                {
                    double p0, q0;
                    BesselAsympt0(x, out p0, out q0);
                    var nn = x - Math.PI / 4;
                    return Math.Sqrt(2 / pi / x) * (p0 * Math.Sin(nn) + q0 * Math.Cos(nn));
                }

                var xsq = x * x;

                var p4 = -41370.35497933148554125235152;
                p4 = 59152134.65686889654273830069 + xsq * p4;
                p4 = -34363712229.79040378171030138 + xsq * p4;
                p4 = 10255208596863.94284509167421 + xsq * p4;
                p4 = -1648605817185729.473122082537 + xsq * p4;
                p4 = 137562431639934407.8571335453 + xsq * p4;
                p4 = -5247065581112764941.297350814 + xsq * p4;
                p4 = 65874732757195549259.99402049 + xsq * p4;
                p4 = -27502866786291095837.01933175 + xsq * p4;

                var q4 = 1.0;
                q4 = 1282.452772478993804176329391 + xsq * q4;
                q4 = 1001702.641288906265666651753 + xsq * q4;
                q4 = 579512264.0700729537480087915 + xsq * q4;
                q4 = 261306575504.1081249568482092 + xsq * q4;
                q4 = 91620380340751.85262489147968 + xsq * q4;
                q4 = 23928830434997818.57439356652 + xsq * q4;
                q4 = 4192417043410839973.904769661 + xsq * q4;
                q4 = 372645883898616588198.9980 + xsq * q4;

                return p4 / q4 + 2 / pi * J0(x) * Math.Log(x);
            }

            /// <summary>Функция Бесселя второго типа, первого порядка</summary>
            /// <param name="x">Аргумент функции Бесселя второго типа, первого порядка</param>
            /// <returns>Значение функции Бесселя второго типа, первого порядка</returns>
            [System.Diagnostics.DebuggerStepThrough]
            public static double Y1(double x)
            {
                if(x > 8)
                {
                    double p0, q0;
                    BesselAsympt1(x, out p0, out q0);
                    var nn = x - 3 * pi / 4;
                    return Math.Sqrt(2 / pi / x) * (p0 * Math.Sin(nn) + q0 * Math.Cos(nn));
                }

                var xsq = x * x;

                var p4 = -2108847.540133123652824139923;
                p4 = 3639488548.124002058278999428 + xsq * p4;
                p4 = -2580681702194.450950541426399 + xsq * p4;
                p4 = 956993023992168.3481121552788 + xsq * p4;
                p4 = -196588746272214065.8820322248 + xsq * p4;
                p4 = 21931073399177975921.11427556 + xsq * p4;
                p4 = -1212297555414509577913.561535 + xsq * p4;
                p4 = 26554738314348543268942.48968 + xsq * p4;
                p4 = -99637534243069222259967.44354 + xsq * p4;

                var q4 = 1.0;
                q4 = 1612.361029677000859332072312 + xsq * q4;
                q4 = 1563282.754899580604737366452 + xsq * q4;
                q4 = 1128686837.169442121732366891 + xsq * q4;
                q4 = 646534088126.5275571961681500 + xsq * q4;
                q4 = 297663212564727.6729292742282 + xsq * q4;
                q4 = 108225825940881955.2553850180 + xsq * q4;
                q4 = 29549879358971486742.90758119 + xsq * q4;
                q4 = 5435310377188854170800.653097 + xsq * q4;
                q4 = 508206736694124324531442.4152 + xsq * q4;

                return x * p4 / q4 + 2 / pi * (J1(x) * Math.Log(x) - 1 / x);
            }

            /// <summary>Функция Бесселя второго типа, n порядка</summary>
            /// <param name="n">Порядок функции Бесселя второго типа</param>
            /// <param name="x">Аргумент функции Бесселя второго типа, n порядка</param>
            /// <returns>Значение функции Бесселя второго типа, n порядка</returns>
            [System.Diagnostics.DebuggerStepThrough]
            public static double Yn(int n, double x)
            {
                double s = 1;

                if(n < 0)
                    if((n = -n) % 2 != 0)
                        s = -1;

                if(n == 0)
                    return Y0(x);

                if(n == 1)
                    return s * Y1(x);

                var a = Y0(x);
                var b = Y1(x);
                for(var i = 1; i <= n - 1; i++)
                {
                    var tmp = b;
                    b = 2 * i / x * b - a;
                    a = tmp;
                }
                return s * b;
            }

            /// <summary>Модифицированная функция Бесселя нулевого порядка </summary>
            /// <param name="x">Аргумент модифицированной функции Бесселя нулевого порядка</param>
            /// <returns>Значение модифицированной функции Бесселя нулевого порядка</returns>
            [System.Diagnostics.DebuggerStepThrough]
            public static double I0(double x)
            {
                if(x < 0) x = -x;

                double b0 = 0;
                double b1 = 0;
                double b2 = 0;

                if(x <= 8)
                {
                    var y = .5 * x - 2;
                    BesselMFirstCheb(-4.41534164647933937950E-18, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 3.33079451882223809783E-17, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -2.43127984654795469359E-16, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.71539128555513303061E-15, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -1.16853328779934516808E-14, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 7.67618549860493561688E-14, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -4.85644678311192946090E-13, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 2.95505266312963983461E-12, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -1.72682629144155570723E-11, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 9.67580903537323691224E-11, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -5.18979560163526290666E-10, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 2.65982372468238665035E-9, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -1.30002500998624804212E-8, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 6.04699502254191894932E-8, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -2.67079385394061173391E-7, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.11738753912010371815E-6, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -4.41673835845875056359E-6, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.64484480707288970893E-5, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -5.75419501008210370398E-5, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.88502885095841655729E-4, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -5.76375574538582365885E-4, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.63947561694133579842E-3, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -4.32430999505057594430E-3, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.05464603945949983183E-2, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -2.37374148058994688156E-2, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 4.93052842396707084878E-2, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -9.49010970480476444210E-2, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.71620901522208775349E-1, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -3.04682672343198398683E-1, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 6.76795274409476084995E-1, ref b0, ref b1, ref b2);
                    return Math.Exp(x) * (.5 * (b0 - b2));
                }

                var z = 32 / x - 2;
                BesselMFirstCheb(-7.23318048787475395456E-18, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -4.83050448594418207126E-18, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 4.46562142029675999901E-17, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 3.46122286769746109310E-17, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -2.82762398051658348494E-16, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -3.42548561967721913462E-16, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 1.77256013305652638360E-15, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 3.81168066935262242075E-15, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -9.55484669882830764870E-15, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -4.15056934728722208663E-14, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 1.54008621752140982691E-14, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 3.85277838274214270114E-13, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 7.18012445138366623367E-13, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -1.79417853150680611778E-12, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -1.32158118404477131188E-11, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -3.14991652796324136454E-11, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 1.18891471078464383424E-11, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 4.94060238822496958910E-10, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 3.39623202570838634515E-9, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 2.26666899049817806459E-8, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 2.04891858946906374183E-7, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 2.89137052083475648297E-6, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 6.88975834691682398426E-5, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 3.36911647825569408990E-3, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 8.04490411014108831608E-1, ref b0, ref b1, ref b2);
                return Math.Exp(x) * (.5 * (b0 - b2)) / Math.Sqrt(x);
            }

            /// <summary>Модифицированная функция Бесселя первого порядка </summary>
            /// <param name="x">Аргумент модифицированной функции Бесселя первого порядка</param>
            /// <returns>Значение модифицированной функции Бесселя первого порядка</returns>
            [System.Diagnostics.DebuggerStepThrough]
            public static double I1(double x)
            {
                double y;
                double b0 = 0;
                double b1 = 0;
                double b2 = 0;

                var z = Math.Abs(x);
                if(z <= 8)
                {
                    y = z / 2 - 2;
                    BesselM1FirstCheb(2.77791411276104639959E-18, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -2.11142121435816608115E-17, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 1.55363195773620046921E-16, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.10559694773538630805E-15, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 7.60068429473540693410E-15, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -5.04218550472791168711E-14, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 3.22379336594557470981E-13, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.98397439776494371520E-12, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 1.17361862988909016308E-11, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -6.66348972350202774223E-11, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 3.62559028155211703701E-10, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.88724975172282928790E-9, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 9.38153738649577178388E-9, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -4.44505912879632808065E-8, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 2.00329475355213526229E-7, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -8.56872026469545474066E-7, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 3.47025130813767847674E-6, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.32731636560394358279E-5, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 4.78156510755005422638E-5, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.61760815825896745588E-4, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 5.12285956168575772895E-4, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.51357245063125314899E-3, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 4.15642294431288815669E-3, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.05640848946261981558E-2, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 2.47264490306265168283E-2, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -5.29459812080949914269E-2, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 1.02643658689847095384E-1, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.76416518357834055153E-1, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 2.52587186443633654823E-1, ref b0, ref b1, ref b2);
                    z = (.5 * (b0 - b2)) * z * Math.Exp(z);
                }
                else
                {
                    y = 32 / z - 2;
                    BesselM1FirstCheb(7.51729631084210481353E-18, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 4.41434832307170791151E-18, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -4.65030536848935832153E-17, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -3.20952592199342395980E-17, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 2.96262899764595013876E-16, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 3.30820231092092828324E-16, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.88035477551078244854E-15, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -3.81440307243700780478E-15, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 1.04202769841288027642E-14, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 4.27244001671195135429E-14, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -2.10154184277266431302E-14, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -4.08355111109219731823E-13, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -7.19855177624590851209E-13, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 2.03562854414708950722E-12, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 1.41258074366137813316E-11, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 3.25260358301548823856E-11, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.89749581235054123450E-11, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -5.58974346219658380687E-10, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -3.83538038596423702205E-9, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -2.63146884688951950684E-8, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -2.51223623787020892529E-7, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -3.88256480887769039346E-6, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.10588938762623716291E-4, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -9.76109749136146840777E-3, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 7.78576235018280120474E-1, ref b0, ref b1, ref b2);
                    z = (.5 * (b0 - b2)) * Math.Exp(z) / Math.Sqrt(z);
                }

                return x < 0 ? -z : z;
            }

            /// <summary>Модифицированная функция Бесселя второго типа, нулевого порядка </summary>
            /// <param name="x">Аргумент модифицированной функции Бесселя второго типа, нулевого порядка</param>
            /// <returns>Значение модифицированной функции Бесселя второго типа, нулевого порядка</returns>
            [System.Diagnostics.DebuggerStepThrough]
            public static double K0(double x)
            {
                Contract.Requires<ArgumentOutOfRangeException>(x > 0, "Функция K0 определена на x > 0");
                //if(x <= 0)
                //    throw new ArgumentOutOfRangeException("x", "Функция K0 определена на x > 0");

                double b0 = 0;
                double b1 = 0;
                double b2 = 0;

                if(x <= 2)
                {
                    var y = x * x - 2;
                    BesselMFirstCheb(1.37446543561352307156E-16, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 4.25981614279661018399E-14, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.03496952576338420167E-11, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.90451637722020886025E-9, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 2.53479107902614945675E-7, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 2.28621210311945178607E-5, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 1.26461541144692592338E-3, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 3.59799365153615016266E-2, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, 3.44289899924628486886E-1, ref b0, ref b1, ref b2);
                    BesselMnExtCheb(y, -5.35327393233902768720E-1, ref b0, ref b1, ref b2);
                    return (.5 * (b0 - b2)) - Math.Log(.5 * x) * I0(x);
                }

                var z = 8 / x - 2;
                BesselMFirstCheb(5.30043377268626276149E-18, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -1.64758043015242134646E-17, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 5.21039150503902756861E-17, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -1.67823109680541210385E-16, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 5.51205597852431940784E-16, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -1.84859337734377901440E-15, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 6.34007647740507060557E-15, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -2.22751332699166985548E-14, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 8.03289077536357521100E-14, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -2.98009692317273043925E-13, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 1.14034058820847496303E-12, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -4.51459788337394416547E-12, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 1.85594911495471785253E-11, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -7.95748924447710747776E-11, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 3.57739728140030116597E-10, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -1.69753450938905987466E-9, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 8.57403401741422608519E-9, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -4.66048989768794782956E-8, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 2.76681363944501510342E-7, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -1.83175552271911948767E-6, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 1.39498137188764993662E-5, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -1.28495495816278026384E-4, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 1.56988388573005337491E-3, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, -3.14481013119645005427E-2, ref b0, ref b1, ref b2);
                BesselMnExtCheb(z, 2.44030308206595545468E0, ref b0, ref b1, ref b2);
                return (.5 * (b0 - b2)) * Math.Exp(-x) / Math.Sqrt(x);
            }

            /// <summary>Модифицированная функция Бесселя второго типа, первого порядка </summary>
            /// <param name="x">Аргумент модифицированной функции Бесселя второго типа, первого порядка</param>
            /// <returns>Значение модифицированной функции Бесселя второго типа, первого порядка</returns>
            /// <exception cref="ArgumentOutOfRangeException">x меньше 0</exception>
            [System.Diagnostics.DebuggerStepThrough]
            public static double K1(double x)
            {
                Contract.Requires<ArgumentOutOfRangeException>(x > 0, "Функция K1 определена на x > 0");
                //if(x <= 0)
                //    throw new ArgumentOutOfRangeException("x", "Функция K1 определена на x > 0");

                double b0 = 0;
                double b1 = 0;
                double b2 = 0;

                if(x <= 2)
                {
                    var y = x * x - 2;
                    BesselM1FirstCheb(-7.02386347938628759343E-18, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -2.42744985051936593393E-15, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -6.66690169419932900609E-13, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.41148839263352776110E-10, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -2.21338763073472585583E-8, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -2.43340614156596823496E-6, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.73028895751305206302E-4, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -6.97572385963986435018E-3, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -1.22611180822657148235E-1, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, -3.53155960776544875667E-1, ref b0, ref b1, ref b2);
                    BesselM1NextCheb(y, 1.52530022733894777053E0, ref b0, ref b1, ref b2);
                    return Math.Log(.5 * x) * I1(x) + (.5 * (b0 - b2)) / x;
                }

                var z = 8 / x - 2;
                BesselM1FirstCheb(-5.75674448366501715755E-18, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 1.79405087314755922667E-17, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -5.68946255844285935196E-17, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 1.83809354436663880070E-16, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -6.05704724837331885336E-16, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 2.03870316562433424052E-15, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -7.01983709041831346144E-15, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 2.47715442448130437068E-14, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -8.97670518232499435011E-14, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 3.34841966607842919884E-13, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -1.28917396095102890680E-12, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 5.13963967348173025100E-12, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -2.12996783842756842877E-11, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 9.21831518760500529508E-11, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -4.19035475934189648750E-10, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 2.01504975519703286596E-9, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -1.03457624656780970260E-8, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 5.74108412545004946722E-8, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -3.50196060308781257119E-7, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 2.40648494783721712015E-6, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -1.93619797416608296024E-5, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 1.95215518471351631108E-4, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, -2.85781685962277938680E-3, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 1.03923736576817238437E-1, ref b0, ref b1, ref b2);
                BesselM1NextCheb(z, 2.72062619048444266945E0, ref b0, ref b1, ref b2);
                return Math.Exp(-x) * (.5 * (b0 - b2)) / Math.Sqrt(x);
            }

            /// <summary>Модифицированная функция Бесселя второго типа, n порядка </summary>
            /// <param name="n">Порядок модифицированной функции Бесселя второго типа</param>
            /// <param name="x">Аргумент модифицированной функции Бесселя второго типа, n порядка</param>
            /// <returns>Значение модифицированной функции Бесселя второго типа, n порядка</returns>
            /// <exception cref="ArgumentOutOfRangeException">при x меньше, либо = 0</exception>
            /// <exception cref="ArgumentOutOfRangeException">|n| больше 31</exception>
            /// <exception cref="OverflowException"></exception>
            [System.Diagnostics.DebuggerStepThrough]
            public static double Kn(int n, double x)
            {
                Contract.Requires<ArgumentOutOfRangeException>(Math.Abs(n) <= 31, "Порядок Kn больше 31");
                Contract.Requires<ArgumentOutOfRangeException>(x > 0, "Область определения x > 0");
                //if(nn > 31)
                //    throw new ArgumentOutOfRangeException("n", "Порядок Kn больше 31"); 
                //if(x <= 0)
                //    throw new ArgumentOutOfRangeException("x", "Область определения x > 0");

                var nn = Math.Abs(n);
                double k;
                double nk1f;
                double t;
                double s;
                double z0;
                double z;
                double fn;
                double pn;
                double pk;
                int i;

                const double eul = 5.772156649015328606065e-1;


                if(x <= 9.55)
                {
                    var ans = .0;
                    z0 = .25 * x * x;
                    fn = 1;
                    pn = 0;
                    var zmn = 1.0;
                    var tox = 2 / x;

                    if(nn > 0)
                    {
                        pn = -eul;
                        k = 1.0;
                        for(i = 1; i < nn; i++)
                        {
                            pn += 1 / k;
                            k++;
                            fn *= k;
                        }

                        zmn = tox;
                        if(nn == 1)
                            ans = 1 / x;
                        else
                        {
                            nk1f = fn / nn;
                            var kf = 1.0;
                            s = nk1f;
                            z = -z0;
                            var zn = 1.0;

                            for(i = 1; i < nn; i++)
                            {
                                nk1f /= nn - i;
                                kf *= i;
                                zn *= z;
                                t = nk1f * zn / kf;
                                s += t;

                                Contract.Assert(c_MaxRealNumber - Math.Abs(t) > Math.Abs(s), "Overflow");
                                if(c_MaxRealNumber - Math.Abs(t) <= Math.Abs(s))
                                    throw new OverflowException();

                                Contract.Assert(tox <= 1 || c_MaxRealNumber / tox >= zmn, "Overflow");
                                if(tox > 1 & c_MaxRealNumber / tox < zmn)
                                    throw new OverflowException();

                                zmn *= tox;
                            }

                            s *= .5;
                            t = Math.Abs(s);

                            Contract.Assert(zmn <= 1 || c_MaxRealNumber / zmn >= t, "Overflow");
                            if(zmn > 1 & c_MaxRealNumber / zmn < t)
                                throw new OverflowException();

                            Contract.Assert(t <= 1 || c_MaxRealNumber / t >= zmn, "Overflow");
                            if(t > 1 && c_MaxRealNumber / t < zmn)
                                throw new OverflowException();

                            ans = s * zmn;
                        }
                    }

                    var tlg = 2 * Math.Log(.5 * x);
                    pk = -eul;

                    if(nn == 0)
                    {
                        pn = pk;
                        t = 1;
                    }
                    else
                    {
                        pn += 1.0 / nn;
                        t = 1 / fn;
                    }

                    s = (pk + pn - tlg) * t;
                    k = 1;

                    do
                    {
                        t *= z0 / (k * (k + nn));
                        pk += 1 / k;
                        pn += 1 / (k + nn);
                        s += (pk + pn - tlg) * t;
                        k++;
                    } while(Math.Abs(t / s) > Eps);

                    s *= .5 / zmn;
                    if(nn % 2 != 0)
                        s = -s;

                    return ans + s;
                }
                if(x > Math.Log(c_MaxRealNumber))
                    return 0;

                k = nn;
                pn = 4 * k * k;
                z0 = 8 * x;
                pk = fn = s = t = 1;
                var nkf = c_MaxRealNumber;
                i = 0;

                do
                {
                    z = pn - pk * pk;
                    t *= z / (fn * z0);
                    nk1f = Math.Abs(t);

                    if(i >= nn & nk1f > nkf)
                        break;

                    nkf = nk1f;
                    s += t;
                    fn++;
                    pk += 2;
                    i++;
                } while(Math.Abs(t / s) > Eps);

                return Math.Exp(-x) * Math.Sqrt(pi / (2 * x)) * s;
            }
        }

        public static class Distribution
        {
            /// <summary>Нормальное распределение</summary>
            [Copyright("1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier")]
            public static class Normal
            {
                /// <summary>Интеграл от exp{-t^2} от нуля до x / .5 sqrt(pi)</summary>
                /// <param name="x">Аргумент функции нормального распределения</param>
                /// <returns>Значение нормального распределения</returns>
                [System.Diagnostics.DebuggerStepThrough]
                public static double ErrorFunction(double x)
                {
                    var s = Math.Sign(x);
                    x = Math.Abs(x);
                    if(x < .5)
                    {
                        var xsq = x * x;
                        var p = .007547728033418631287834;
                        p = .288805137207594084924010 + xsq * p;
                        p = 14.3383842191748205576712 + xsq * p;
                        p = 38.0140318123903008244444 + xsq * p;
                        p = 3017.82788536507577809226 + xsq * p;
                        p = 7404.07142710151470082064 + xsq * p;
                        p = 80437.3630960840172832162 + xsq * p;

                        var q = .0;
                        q = 1.00000000000000000000000 + xsq * q;
                        q = 38.0190713951939403753468 + xsq * q;
                        q = 658.070155459240506326937 + xsq * q;
                        q = 6379.60017324428279487120 + xsq * q;
                        q = 34216.5257924628539769006 + xsq * q;
                        q = 80437.3630960840172826266 + xsq * q;
                        return s * 1.1283791670955125738961589031 * x * p / q;
                    }
                    return x >= 10 ? s : s * (1 - ErrorFunctionComform(x));
                }

                [System.Diagnostics.DebuggerStepThrough]
                public static double ErrorFunctionComform(double x)
                {
                    if(x < 0)
                        return 2 - ErrorFunctionComform(-x);

                    if(x < .5)
                        return 1 - ErrorFunction(x);

                    if(x >= 10)
                        return 0;

                    var p = .0;
                    p = .5641877825507397413087057563 + x * p;
                    p = 9.675807882987265400604202961 + x * p;
                    p = 77.08161730368428609781633646 + x * p;
                    p = 368.5196154710010637133875746 + x * p;
                    p = 1143.262070703886173606073338 + x * p;
                    p = 2320.439590251635247384768711 + x * p;
                    p = 2898.0293292167655611275846 + x * p;
                    p = 1826.3348842295112592168999 + x * p;

                    var q = 1.0;
                    q = 17.14980943627607849376131193 + x * q;
                    q = 137.1255960500622202878443578 + x * q;
                    q = 661.7361207107653469211984771 + x * q;
                    q = 2094.384367789539593790281779 + x * q;
                    q = 4429.612803883682726711528526 + x * q;
                    q = 6089.5424232724435504633068 + x * q;
                    q = 4958.82756472114071495438422 + x * q;
                    q = 1826.3348842295112595576438 + x * q;
                    return Math.Exp(-(x * x)) * p / q;
                }

                [System.Diagnostics.DebuggerStepThrough]
                public static double NormalDistribution(double x)
                {
                    return .5 * (ErrorFunction(x / 1.41421356237309504880) + 1);
                }

                [System.Diagnostics.DebuggerStepThrough]
                public static double ErrorFunctionInversed(double e)
                {
                    return NormalDistributionInversed(.5 * (e + 1)) / sqrt2;
                }

                [System.Diagnostics.DebuggerStepThrough]
                public static double NormalDistributionInversed(double y0)
                {
                    const double expm2 = .13533528323661269189;
                    const double lc_S2Pi = 2.50662827463100050242;

                    if(y0 <= 0)
                        return -c_MaxRealNumber;

                    if(y0 >= 1)
                        return c_MaxRealNumber;

                    var code = 1;
                    var y = y0;
                    if(y > 1 - expm2)
                    {
                        y = 1 - y;
                        code = 0;
                    }

                    if(y > expm2)
                    {
                        y -= .5;
                        var y2 = y * y;
                        var p0 = -59.9633501014107895267;
                        p0 = 98.0010754185999661536 + y2 * p0;
                        p0 = -56.6762857469070293439 + y2 * p0;
                        p0 = 13.9312609387279679503 + y2 * p0;
                        p0 = -1.23916583867381258016 + y2 * p0;

                        var q0 = 1d;
                        q0 = 1.95448858338141759834 + y2 * q0;
                        q0 = 4.67627912898881538453 + y2 * q0;
                        q0 = 86.3602421390890590575 + y2 * q0;
                        q0 = -225.462687854119370527 + y2 * q0;
                        q0 = 200.260212380060660359 + y2 * q0;
                        q0 = -82.0372256168333339912 + y2 * q0;
                        q0 = 15.9056225126211695515 + y2 * q0;
                        q0 = -1.18331621121330003142 + y2 * q0;
                        return (y + y * y2 * p0 / q0) * lc_S2Pi;
                    }

                    var x = Math.Sqrt(-2 * Math.Log(y));
                    var x0 = x - Math.Log(x) / x;
                    var z = 1 / x;

                    double x1;
                    if(x < 8)
                    {
                        var p1 = 4.05544892305962419923;
                        p1 = 31.5251094599893866154 + z * p1;
                        p1 = 57.1628192246421288162 + z * p1;
                        p1 = 44.0805073893200834700 + z * p1;
                        p1 = 14.6849561928858024014 + z * p1;
                        p1 = 2.18663306850790267539 + z * p1;
                        p1 = -(1.40256079171354495875 * .1) + z * p1;
                        p1 = -(3.50424626827848203418 * .01) + z * p1;
                        p1 = -(8.57456785154685413611 * .0001) + z * p1;

                        var q1 = 1.0;
                        q1 = 15.7799883256466749731 + z * q1;
                        q1 = 45.3907635128879210584 + z * q1;
                        q1 = 41.3172038254672030440 + z * q1;
                        q1 = 15.0425385692907503408 + z * q1;
                        q1 = 2.50464946208309415979 + z * q1;
                        q1 = -(1.42182922854787788574 * .1) + z * q1;
                        q1 = -(3.80806407691578277194 * .01) + z * q1;
                        q1 = -(9.33259480895457427372 * .0001) + z * q1;
                        x1 = z * p1 / q1;
                    }
                    else
                    {
                        var p2 = 3.23774891776946035970;
                        p2 = 6.91522889068984211695 + z * p2;
                        p2 = 3.93881025292474443415 + z * p2;
                        p2 = 1.33303460815807542389 + z * p2;
                        p2 = 2.01485389549179081538 * .1 + z * p2;
                        p2 = 1.23716634817820021358 * .01 + z * p2;
                        p2 = 3.01581553508235416007 * .0001 + z * p2;
                        p2 = 2.65806974686737550832 * .000001 + z * p2;
                        p2 = 6.23974539184983293730 * .000000001 + z * p2;

                        var q2 = 1d;
                        q2 = 6.02427039364742014255 + z * q2;
                        q2 = 3.67983563856160859403 + z * q2;
                        q2 = 1.37702099489081330271 + z * q2;
                        q2 = 2.16236993594496635890 * .1 + z * q2;
                        q2 = 1.34204006088543189037 * .01 + z * q2;
                        q2 = 3.28014464682127739104 * .0001 + z * q2;
                        q2 = 2.89247864745380683936 * .000001 + z * q2;
                        q2 = 6.79019408009981274425 * .000000001 + z * q2;
                        x1 = z * p2 / q2;
                    }
                    return code == 0 ? x0 - x1 : x1 - x0;
                }
            }

            public static class Student
            {
                [System.Diagnostics.DebuggerStepThrough]
                public static double StudenttDistribution(int k, double t)
                {
                    Contract.Requires<ArgumentOutOfRangeException>(k > 0, "Функция определена для положительных чисел k");
                    //if(k <= 0)
                    //    throw new ArgumentOutOfRangeException("k", "Функция определна для положительных чисел k");

                    if(Math.Abs(t - 0) < Eps) return .5;
                    if(t < -2)
                        return .5 * IncompliteBetaFunction.IncompleteBeta(.5 * k, .5, k / (k + t * t));

                    var x = t < 0 ? -t : t;

                    double rk = k;
                    var z = 1 + x * x / rk;
                    double tz;
                    double f;
                    double p;
                    int j;
                    if(k % 2 != 0)
                    {
                        var xsqk = x / Math.Sqrt(rk);
                        p = Math.Atan(xsqk);
                        if(k > 1)
                        {
                            f = 1;
                            tz = 1;
                            j = 3;
                            while(j <= k - 2 & tz / f > Eps)
                            {
                                tz *= (j - 1) / (z * j);
                                f += tz;
                                j += 2;
                            }
                            p += f * xsqk / z;
                        }
                        p *= 2 / pi;
                    }
                    else
                    {
                        f = tz = 1;
                        j = 2;
                        while(j <= k - 2 & tz / f > Eps)
                        {
                            tz *= (j - 1) / (z * j);
                            f += tz;
                            j += 2;
                        }
                        p = f * x / Math.Sqrt(z * rk);
                    }

                    return .5 + .5 * (t < 0 ? -p : p);
                }

                [System.Diagnostics.DebuggerStepThrough]
                public static double StudenttDistributionInversed(int k, double p)
                {
                    Contract.Requires<ArgumentOutOfRangeException>(k > 0, "k должно быть больше 0");
                    Contract.Requires<ArgumentOutOfRangeException>(p > 0 && p < 1, "p должно быть больше 0 и меньше 1");
                    //if(k <= 0)
                    //    throw new ArgumentOutOfRangeException("k", "k должно быть больше 0");
                    //if(p <= 0 || p >= 1)
                    //    throw new ArgumentOutOfRangeException("p", "p должно быть больше 0 и меньше единицы");

                    double z;

                    double rk = k;
                    if(p > .25 && p < .75)
                    {
                        if(Math.Abs(p - .5) < Eps) return 0;
                        z = IncompliteBetaFunction.IncompleteBetaInversed(.5, .5 * rk, Math.Abs(1 - 2 * p));
                        var t = Math.Sqrt(rk * z / (1 - z));
                        return p < 0 ? -t : t;
                    }

                    var rflg = -1;
                    if(p >= .5)
                    {
                        p = 1 - p;
                        rflg = 1;
                    }

                    z = IncompliteBetaFunction.IncompleteBetaInversed(.5 * rk, .5, 2 * p);
                    return c_MaxRealNumber * z < rk ? rflg * c_MaxRealNumber : rflg * Math.Sqrt(rk / z - rk);
                }

                [System.Diagnostics.DebuggerStepThrough]
                public static double QuantileHi2(int n, double alpha = .05)
                {
                    if(alpha < .001 || alpha > .999)
                        throw new NotSupportedException("Значения alpha < 0.001 и > 0.999 не поддерживаются");

                    var sqrt_n = Math.Sqrt(n);
                    var b = (alpha < .5)
                                ? -2.0637 * Math.Pow(Math.Log(1 / alpha) - .16, .4274) + 1.5774
                                : 2.0637 * Math.Pow(Math.Log(1 / (1 - alpha)) - .16, .4274) - 1.5774;

                    var A = b / sqrt2;
                    var B = 2 * (b * b - 1) / 3;
                    var C = b * (b * b - 7) / (9 * sqrt2);
                    var D = ((6 * b * b * b + 14 * b) * b - 32) / 405;
                    var E = b * ((9 * b * b * b + 256 * b) * b - 433) / (4826 * sqrt2);
                    return n + A * sqrt_n + B + C / sqrt_n + D / n + E / n / sqrt_n;
                }

            }
        }

        public static class IncompliteBetaFunction
        {
            /*************************************************************************
            Incomplete beta integral

            Returns incomplete beta integral of the arguments, evaluated
            from zero to x.  The function is defined as

                             x
                -            -
               | (a+b)      | |  a-1     b-1
             -----------    |   t   (1-t)   dt.
              -     -     | |
             | (a) | (b)   -
                            0

            The domain of definition is 0 <= x <= 1.  In this
            implementation a and b are restricted to positive values.
            The integral from x to 1 may be obtained by the symmetry
            relation

               1 - incbet( a, b, x )  =  incbet( b, a, 1-x ).

            The integral is evaluated by a continued fraction expansion
            or, when b*x is small, by a power series.

            ACCURACY:

            Tested at uniformly distributed random points (a,b,x) with a and b
            in "domain" and x between 0 and 1.
                                                   Relative error
            arithmetic   domain     # trials      peak         rms
               IEEE      0,5         10000       6.9e-15     4.5e-16
               IEEE      0,85       250000       2.2e-13     1.7e-14
               IEEE      0,1000      30000       5.3e-12     6.3e-13
               IEEE      0,10000    250000       9.3e-11     7.1e-12
               IEEE      0,100000    10000       8.7e-10     4.8e-11
            Outputs smaller than the IEEE gradual underflow threshold
            were excluded from these statistics.

            Cephes Math Library, Release 2.8:  June, 2000
            Copyright 1984, 1995, 2000 by Stephen L. Moshier
            *************************************************************************/
            [System.Diagnostics.DebuggerStepThrough]
            public static double IncompleteBeta(double a, double b, double x)
            {
                Contract.Assert(a > 0, "a > 0");
                if(a <= 0)
                    throw new ArgumentOutOfRangeException("a", "a должна быть > 0");

                Contract.Assert(b > 0, "b > 0");
                if(b <= 0)
                    throw new ArgumentOutOfRangeException("b", "b должна быть > 0");

                Contract.Assert(x >= 0 && x <= 1, "0 <= x <= 1");
                if(x < 0 || x > 1)
                    throw new ArgumentOutOfRangeException("x", "a должна быть в пределах [0;1]");

                const double big = 4.503599627370496e15;
                const double biginv = 2.22044604925031308085e-16;
                const double maxgam = 171.624376956302725;
                const double minlog = c_LogMinRealNumber;
                const double maxlog = c_LogMaxRealNumber;

                if(Math.Abs(x - 0) < Eps) return 0;
                if(Math.Abs(x - 1) < Eps) return 1;
                var flag = 0;

                if(b * x <= 1 && x <= .95)
                    return IncompleteBetaPowerSeries(a, b, x, maxgam);

                var w = 1 - x;

                double t;
                double xc;
                if(x <= a / (a + b)) xc = w;
                else
                {
                    flag = 1;
                    t = a;
                    a = b;
                    b = t;
                    xc = x;
                    x = w;
                }

                if(flag == 1 && b * x <= 1 && x <= .95)
                {
                    t = IncompleteBetaPowerSeries(a, b, x, maxgam);
                    return 1 - (t <= Eps ? Eps : t);
                }

                var y = x * (a + b - 2) - (a - 1);
                w = y < 0 ? IncompleteBetaFractionExpansion(a, b, x, big, biginv) : IncompleteBetaFractionExpansion2(a, b, x, big, biginv) / xc;

                y = a * Math.Log(x);
                t = b * Math.Log(xc);

                if(a + b < maxgam && Math.Abs(y) < maxlog && Math.Abs(t) < maxlog)
                {
                    t = Math.Pow(xc, b) * Math.Pow(x, a) / a * w * Gamma.G(a + b) / (Gamma.G(a) * Gamma.G(b));
                    return flag == 1 ? 1 - (t <= Eps ? Eps : t) : t;
                }
                int sign;
                y += t + Gamma.LnG(a + b, out sign) - Gamma.LnG(a, out sign) - Gamma.LnG(b, out sign);
                y += Math.Log(w / a);
                t = y < minlog ? 0 : Math.Exp(y);
                if(flag == 1)
                    t = 1 - (t <= Eps ? Eps : t);
                return t;
            }


            /*************************************************************************
            Inverse of imcomplete beta integral

            Given y, the function finds x such that

             incbet( a, b, x ) = y .

            The routine performs interval halving or Newton iterations to find the
            root of incbet(a,b,x) - y = 0.


            ACCURACY:

                                 Relative error:
                           x     a,b
            arithmetic   domain  domain  # trials    peak       rms
               IEEE      0,1    .5,10000   50000    5.8e-12   1.3e-13
               IEEE      0,1   .25,100    100000    1.8e-13   3.9e-15
               IEEE      0,1     0,5       50000    1.1e-12   5.5e-15
            With a and b constrained to half-integer or integer values:
               IEEE      0,1    .5,10000   50000    5.8e-12   1.1e-13
               IEEE      0,1    .5,100    100000    1.7e-14   7.9e-16
            With a = .5, b constrained to half-integer or integer values:
               IEEE      0,1    .5,10000   10000    8.3e-11   1.0e-11

            Cephes Math Library Release 2.8:  June, 2000
            Copyright 1984, 1996, 2000 by Stephen L. Moshier
            *************************************************************************/

            [System.Diagnostics.DebuggerStepThrough]
            public static double IncompleteBetaInversed(double a, double b, double y)
            {
                Contract.Requires<ArgumentOutOfRangeException>(y >= 0 && y <= 1, "y должен находиться в пределах [0;1]");
                //if(y < 0 || y > 1)
                //    throw new ArgumentOutOfRangeException("y", "y должен находиться в пределах [0;1]");

                //
                // special cases
                //
                if(Math.Abs(y) < Eps) return 0;
                if(Math.Abs(y - 1) < Eps) return 1;

                //
                // these initializations are not really necessary,
                // but without them compiler complains about 'possibly uninitialized variables'.
                //
                double dithresh = 0;
                var rflg = 0;
                double aaa = 0;
                double bbb = 0;
                double y0 = 0;
                double x = 0;
                double yyy = 0;
                double lgm = 0;
                var dir = 0;
                double di = 0;

                //
                // normal initializations
                //
                double x0 = 0;
                double yl = 0;
                double x1 = 1;
                double yh = 1;
                var nflg = 0;
                var mainlooppos = 0;
                const int ihalve = 1;
                const int ihalvecycle = 2;
                const int newt = 3;
                const int newtcycle = 4;
                const int breaknewtcycle = 5;
                const int breakihalvecycle = 6;

                //
                // main loop
                //
                while(true)
                {
                    //
                    // start
                    //
                    double d;
                    double yp;
                    if(mainlooppos == 0)
                    {
                        if(a <= 1 || b <= 1)
                        {
                            dithresh = 1e-6;
                            rflg = 0;
                            aaa = a;
                            bbb = b;
                            y0 = y;
                            x = aaa / (aaa + bbb);
                            yyy = IncompleteBeta(aaa, bbb, x);
                            mainlooppos = ihalve;
                            continue;
                        }

                        dithresh = 1e-4;
                        yp = -Distribution.Normal.NormalDistributionInversed(y);

                        if(y > .5)
                        {
                            rflg = 1;
                            aaa = b;
                            bbb = a;
                            y0 = 1 - y;
                            yp = -yp;
                        }
                        else
                        {
                            rflg = 0;
                            aaa = a;
                            bbb = b;
                            y0 = y;
                        }

                        lgm = (yp * yp - 3) / 6;
                        x = 2 / (1 / (2 * aaa - 1) + 1 / (2 * bbb - 1));
                        d = yp * Math.Sqrt(x + lgm) / x - (1 / (2 * bbb - 1) - 1 / (2 * aaa - 1)) *
                            (lgm + 5 / 6 - 2 / (3 * x));

                        d *= 2;

                        if(d < Math.Log(c_MinRealNumber))
                        {
                            x = 0;
                            break;
                        }

                        x = aaa / (aaa + bbb * Math.Exp(d));
                        yyy = IncompleteBeta(aaa, bbb, x);
                        yp = (yyy - y0) / y0;

                        if(Math.Abs(yp) < .2)
                        {
                            mainlooppos = newt;
                            continue;
                        }
                        mainlooppos = ihalve;
                        continue;
                    }

                    //
                    // ihalve
                    //
                    var i = 0;
                    if(mainlooppos == ihalve)
                    {
                        dir = 0;
                        di = .5;
                        mainlooppos = ihalvecycle;
                        continue;
                    }

                    //
                    // ihalvecycle
                    //
                    if(mainlooppos == ihalvecycle)
                    {
                        if(i <= 99)
                        {
                            if(i != 0)
                            {
                                x = x0 + di * (x1 - x0);

                                if(Math.Abs(x - 1) < Eps)
                                    x = 1 - Eps;

                                if(Math.Abs(x) < Eps)
                                {
                                    di = .5;
                                    x = x0 + di * (x1 - x0);

                                    if(Math.Abs(x) < Eps) break;
                                }

                                yyy = IncompleteBeta(aaa, bbb, x);
                                yp = (x1 - x0) / (x1 + x0);

                                if(Math.Abs(yp) < dithresh)
                                {
                                    mainlooppos = newt;
                                    continue;
                                }

                                yp = (yyy - y0) / y0;

                                if(Math.Abs(yp) < dithresh)
                                {
                                    mainlooppos = newt;
                                    continue;
                                }
                            }

                            if(yyy < y0)
                            {
                                x0 = x;
                                yl = yyy;
                                if(dir >= 0)
                                    di = dir > 3
                                        ? 1 - (1 - di) * (1 - di)
                                        : (dir > 1 ? .5 * di + .5 : (y0 - yyy) / (yh - yl));
                                else
                                {
                                    dir = 0;
                                    di = .5;
                                }
                                dir++;
                                if(x0 > .75)
                                {
                                    if(rflg == 1)
                                    {
                                        rflg = 0;
                                        aaa = a;
                                        bbb = b;
                                        y0 = y;
                                    }
                                    else
                                    {
                                        rflg = 1;
                                        aaa = b;
                                        bbb = a;
                                        y0 = 1 - y;
                                    }
                                    x = 1 - x;
                                    yyy = IncompleteBeta(aaa, bbb, x);
                                    x0 = 0;
                                    yl = 0;
                                    x1 = 1;
                                    yh = 1;
                                    mainlooppos = ihalve;
                                    continue;
                                }
                            }
                            else
                            {
                                x1 = x;
                                if(rflg == 1 && x1 < Eps)
                                {
                                    x = 0;
                                    break;
                                }
                                yh = yyy;

                                if(dir <= 0) di = dir < -3 ? di * di : (dir < -1 ? di * .5 : (yyy - y0) / (yh - yl));
                                else
                                {
                                    dir = 0;
                                    di = .5;
                                }

                                dir--;
                            }
                            i++;
                            mainlooppos = ihalvecycle;
                            continue;
                        }
                        mainlooppos = breakihalvecycle;
                        continue;
                    }

                    //
                    // breakihalvecycle
                    //
                    if(mainlooppos == breakihalvecycle)
                    {
                        if(x0 >= 1)
                        {
                            x = 1 - Eps;
                            break;
                        }
                        if(x <= 0)
                        {
                            x = 0;
                            break;
                        }
                        mainlooppos = newt;
                        continue;
                    }

                    //
                    // newt
                    //
                    if(mainlooppos == newt)
                    {
                        if(nflg != 0) break;

                        nflg = 1;
                        int s;
                        lgm = Gamma.LnG(aaa + bbb, out s) - Gamma.LnG(aaa, out s) - Gamma.LnG(bbb, out s);
                        mainlooppos = newtcycle;
                        continue;
                    }

                    //
                    // newtcycle
                    //
                    if(mainlooppos == newtcycle)
                    {
                        if(i <= 7)
                        {
                            if(i != 0) yyy = IncompleteBeta(aaa, bbb, x);

                            if(yyy < yl)
                            {
                                x = x0;
                                yyy = yl;
                            }
                            else
                            {
                                if(yyy > yh)
                                {
                                    x = x1;
                                    yyy = yh;
                                }
                                else
                                {
                                    if(yyy < y0)
                                    {
                                        x0 = x;
                                        yl = yyy;
                                    }
                                    else
                                    {
                                        x1 = x;
                                        yh = yyy;
                                    }
                                }
                            }
                            if(Math.Abs(x - 1) < Eps || Math.Abs(x) < Eps)
                            {
                                mainlooppos = breaknewtcycle;
                                continue;
                            }

                            d = (aaa - 1) * Math.Log(x) + (bbb - 1) * Math.Log(1 - x) + lgm;

                            if(d < Math.Log(c_MinRealNumber)) break;

                            if(d > Math.Log(c_MaxRealNumber))
                            {
                                mainlooppos = breaknewtcycle;
                                continue;
                            }
                            d = (yyy - y0) / Math.Exp(d);
                            var xt = x - d;

                            if(xt <= x0)
                            {
                                yyy = (x - x0) / (x1 - x0);
                                xt = x0 + .5 * yyy * (x - x0);

                                if(xt <= 0)
                                {
                                    mainlooppos = breaknewtcycle;
                                    continue;
                                }
                            }

                            if(xt >= x1)
                            {
                                yyy = (x1 - x) / (x1 - x0);
                                xt = x1 - .5 * yyy * (x1 - x);
                                if(xt >= 1)
                                {
                                    mainlooppos = breaknewtcycle;
                                    continue;
                                }
                            }

                            x = xt;
                            if(Math.Abs(d / x) < 128 * Eps) break;
                            i++;
                            mainlooppos = newtcycle;
                            continue;
                        }
                        mainlooppos = breaknewtcycle;
                        continue;
                    }

                    //
                    // breaknewtcycle
                    //
                    if(mainlooppos != breaknewtcycle) continue;

                    dithresh = 256 * Eps;
                    mainlooppos = ihalve;
                }

                //
                // done
                //
                return rflg == 0 ? x : (x <= Eps ? 1 - Eps : 1 - x);
            }


            /*************************************************************************
            Continued fraction expansion #1 for incomplete beta integral

            Cephes Math Library, Release 2.8:  June, 2000
            Copyright 1984, 1995, 2000 by Stephen L. Moshier
            *************************************************************************/

            private static double IncompleteBetaFractionExpansion(double a, double b, double x, double big, double biginv)
            {
                var k1 = a;
                var k2 = a + b;
                var k3 = a;
                var k4 = a + 1;
                var k5 = 1.0;
                var k6 = b - 1;
                var k7 = k4;
                var k8 = a + 2;
                var pkm2 = 0.0;
                var qkm2 = 1.0;
                var pkm1 = 1.0;
                var qkm1 = 1.0;
                var ans = 1.0;
                var r = 1.0;
                var n = 0;

                const double thresh = 3 * Eps;
                do
                {
                    var xk = -(x * k1 * k2 / (k3 * k4));
                    var pk = pkm1 + pkm2 * xk;
                    var qk = qkm1 + qkm2 * xk;
                    pkm2 = pkm1;
                    pkm1 = pk;
                    qkm2 = qkm1;
                    qkm1 = qk;
                    xk = x * k5 * k6 / (k7 * k8);
                    pk = pkm1 + pkm2 * xk;
                    qk = qkm1 + qkm2 * xk;
                    pkm2 = pkm1;
                    pkm1 = pk;
                    qkm2 = qkm1;
                    qkm1 = qk;

                    if(Math.Abs(qk) > Eps)
                        r = pk / qk;
                    double t;
                    if(Math.Abs(r) < Eps) t = 1;
                    else
                    {
                        t = Math.Abs((ans - r) / r);
                        ans = r;
                    }

                    if(t < thresh) break;

                    k1 += 1;
                    k2 += 1;
                    k3 += 2;
                    k4 += 2;
                    k5 += 1;
                    k6 += 1;
                    k7 += 2;
                    k8 += 2;

                    if(Math.Abs(qk) + Math.Abs(pk) > big)
                    {
                        pkm2 *= biginv;
                        pkm1 *= biginv;
                        qkm2 *= biginv;
                        qkm1 *= biginv;
                    }

                    if(Math.Abs(qk) >= biginv && Math.Abs(pk) >= biginv) continue;

                    pkm2 *= big;
                    pkm1 *= big;
                    qkm2 *= big;
                    qkm1 *= big;
                } while(++n != 300);
                return ans;
            }


            /*************************************************************************
            Continued fraction expansion #2
            for incomplete beta integral

            Cephes Math Library, Release 2.8:  June, 2000
            Copyright 1984, 1995, 2000 by Stephen L. Moshier
            *************************************************************************/
            private static double IncompleteBetaFractionExpansion2(double a, double b, double x, double big, double biginv)
            {
                var k1 = a;
                var k2 = b - 1;
                var k3 = a;
                var k4 = a + 1;
                var k5 = 1.0;
                var k6 = a + b;
                var k7 = a + 2;
                var k8 = a + 2;
                var pkm2 = 0.0;
                var qkm2 = 1.0;
                var pkm1 = 1.0;
                var qkm1 = 1.0;
                var z = x / (1 - x);
                var ans = 1.0;
                var r = 1.0;
                var n = 0;
                const double thresh = 3 * Eps;

                do
                {
                    var xk = -(z * k1 * k2 / (k3 * k4));
                    var pk = pkm1 + pkm2 * xk;
                    var qk = qkm1 + qkm2 * xk;
                    pkm2 = pkm1;
                    pkm1 = pk;
                    qkm2 = qkm1;
                    qkm1 = qk;
                    xk = z * k5 * k6 / (k7 * k8);
                    pk = pkm1 + pkm2 * xk;
                    qk = qkm1 + qkm2 * xk;
                    pkm2 = pkm1;
                    pkm1 = pk;
                    qkm2 = qkm1;
                    qkm1 = qk;

                    if(Math.Abs(qk) > Eps) r = pk / qk;
                    double t;

                    if(Math.Abs(r) < Eps) t = 1;
                    else
                    {
                        t = Math.Abs((ans - r) / r);
                        ans = r;
                    }

                    if(t < thresh) break;

                    k1 += 1;
                    k2 += 1;
                    k3 += 2;
                    k4 += 2;
                    k5 += 1;
                    k6 += 1;
                    k7 += 2;
                    k8 += 2;

                    if(Math.Abs(qk) + Math.Abs(pk) > big)
                    {
                        pkm2 *= biginv;
                        pkm1 *= biginv;
                        qkm2 *= biginv;
                        qkm1 *= biginv;
                    }
                    if(Math.Abs(qk) >= biginv && Math.Abs(pk) >= biginv) continue;

                    pkm2 *= big;
                    pkm1 *= big;
                    qkm2 *= big;
                    qkm1 *= big;
                } while(++n != 300);
                return ans;
            }


            /*************************************************************************
            Power series for incomplete beta integral.
            Use when b*x is small and x not too close to 1.

            Cephes Math Library, Release 2.8:  June, 2000
            Copyright 1984, 1995, 2000 by Stephen L. Moshier
            *************************************************************************/
            private static double IncompleteBetaPowerSeries(double a, double b, double x, double maxgam)
            {
                var ai = 1 / a;
                var u = (1 - b) * x;
                var v = u / (a + 1);
                var t1 = v;
                var t = u;
                var n = 2.0;
                var s = 0.0;
                var z = Eps * ai;
                while(Math.Abs(v) > z)
                {
                    u = (n - b) * x / n;
                    t *= u;
                    v = t / (a + n);
                    s += v;
                    n++;
                }
                s += t1;
                s += ai;
                u = a * Math.Log(x);
                if((a + b) < maxgam && Math.Abs(u) < Math.Log(c_MaxRealNumber))
                {
                    t = Gamma.G(a + b) / (Gamma.G(a) * Gamma.G(b));
                    s *= t * Math.Pow(x, a);
                }
                else
                {
                    int sign;
                    t = Gamma.LnG(a + b, out sign) - Gamma.LnG(a, out sign) - Gamma.LnG(b, out sign) + u + Math.Log(s);
                    s = t < Math.Log(c_MinRealNumber) ? 0 : Math.Exp(t);
                }
                return s;
            }
        }

        public class TrigonometryFunctionsIntegrals
        {
            /*************************************************************************
            Sine and cosine integrals

            Evaluates the integrals

                                     x
                                     -
                                    |  cos t - 1
              Ci(x) = eul + ln x +  |  --------- dt,
                                    |      t
                                   -
                                    0
                        x
                        -
                       |  sin t
              Si(x) =  |  ----- dt
                       |    t
                      -
                       0

            where eul = 0.57721566490153286061 is Euler's constant.
            The integrals are approximated by rational functions.
            For x > 8 auxiliary functions f(x) and g(x) are employed
            such that

            Ci(x) = f(x) sin(x) - g(x) cos(x)
            Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)


            ACCURACY:
               Test interval = [0,50].
            Absolute error, except relative when > 1:
            arithmetic   function   # trials      peak         rms
               IEEE        Si        30000       4.4e-16     7.3e-17
               IEEE        Ci        30000       6.9e-16     5.1e-17

            Cephes Math Library Release 2.1:  January, 1989
            Copyright 1984, 1987, 1989 by Stephen L. Moshier
            *************************************************************************/

            [System.Diagnostics.DebuggerStepThrough]
            public static void SineCosineIntegrals(double x, out double si, out double ci)
            {

                int sg;
                if(x >= 0) sg = 0;
                else
                {
                    sg = -1;
                    x = -x;
                }

                if(Math.Abs(x) < Eps)
                {
                    si = 0;
                    ci = -c_MaxRealNumber;
                    return;
                }

                if(x > 1E9)
                {
                    si = pi_2 - Math.Cos(x) / x;
                    ci = Math.Sin(x) / x;
                    return;
                }

                double s;
                double c;
                double z;
                if(x <= 4)
                {
                    z = x * x;
                    var sn = -8.39167827910303881427E-11;
                    sn = sn * z + 4.62591714427012837309E-8;
                    sn = sn * z - 9.75759303843632795789E-6;
                    sn = sn * z + 9.76945438170435310816E-4;
                    sn = sn * z - 4.13470316229406538752E-2;
                    sn = sn * z + 1.00000000000000000302E0;

                    var sd = 2.03269266195951942049E-12;
                    sd = sd * z + 1.27997891179943299903E-9;
                    sd = sd * z + 4.41827842801218905784E-7;
                    sd = sd * z + 9.96412122043875552487E-5;
                    sd = sd * z + 1.42085239326149893930E-2;
                    sd = sd * z + 9.99999999999999996984E-1;
                    s = x * sn / sd;

                    var cn = 2.02524002389102268789E-11;
                    cn = cn * z - 1.35249504915790756375E-8;
                    cn = cn * z + 3.59325051419993077021E-6;
                    cn = cn * z - 4.74007206873407909465E-4;
                    cn = cn * z + 2.89159652607555242092E-2;
                    cn = cn * z - 1.00000000000000000080E0;

                    var cd = 4.07746040061880559506E-12;
                    cd = cd * z + 3.06780997581887812692E-9;
                    cd = cd * z + 1.23210355685883423679E-6;
                    cd = cd * z + 3.17442024775032769882E-4;
                    cd = cd * z + 5.10028056236446052392E-2;
                    cd = cd * z + 4.00000000000000000080E0;
                    c = z * cn / cd;
                    si = sg == 0 ? s : -s;
                    ci = 0.57721566490153286061 + Math.Log(x) + c;
                    return;
                }

                s = Math.Sin(x);
                c = Math.Cos(x);
                z = 1 / (x * x);

                double gd;
                double gn;
                double fd;
                double fn;
                double g;
                double f;
                if(x < 8)
                {
                    fn = 4.23612862892216586994E0;
                    fn = fn * z + 5.45937717161812843388E0;
                    fn = fn * z + 1.62083287701538329132E0;
                    fn = fn * z + 1.67006611831323023771E-1;
                    fn = fn * z + 6.81020132472518137426E-3;
                    fn = fn * z + 1.08936580650328664411E-4;
                    fn = fn * z + 5.48900223421373614008E-7;

                    fd = 1.00000000000000000000E0;
                    fd = fd * z + 8.16496634205391016773E0;
                    fd = fd * z + 7.30828822505564552187E0;
                    fd = fd * z + 1.86792257950184183883E0;
                    fd = fd * z + 1.78792052963149907262E-1;
                    fd = fd * z + 7.01710668322789753610E-3;
                    fd = fd * z + 1.10034357153915731354E-4;
                    fd = fd * z + 5.48900252756255700982E-7;
                    f = fn / (x * fd);

                    gn = 8.71001698973114191777E-2;
                    gn = gn * z + 6.11379109952219284151E-1;
                    gn = gn * z + 3.97180296392337498885E-1;
                    gn = gn * z + 7.48527737628469092119E-2;
                    gn = gn * z + 5.38868681462177273157E-3;
                    gn = gn * z + 1.61999794598934024525E-4;
                    gn = gn * z + 1.97963874140963632189E-6;
                    gn = gn * z + 7.82579040744090311069E-9;

                    gd = 1.00000000000000000000E0;
                    gd = gd * z + 1.64402202413355338886E0;
                    gd = gd * z + 6.66296701268987968381E-1;
                    gd = gd * z + 9.88771761277688796203E-2;
                    gd = gd * z + 6.22396345441768420760E-3;
                    gd = gd * z + 1.73221081474177119497E-4;
                    gd = gd * z + 2.02659182086343991969E-6;
                    gd = gd * z + 7.82579218933534490868E-9;
                    g = z * gn / gd;
                }
                else
                {
                    fn = 4.55880873470465315206E-1;
                    fn = fn * z + 7.13715274100146711374E-1;
                    fn = fn * z + 1.60300158222319456320E-1;
                    fn = fn * z + 1.16064229408124407915E-2;
                    fn = fn * z + 3.49556442447859055605E-4;
                    fn = fn * z + 4.86215430826454749482E-6;
                    fn = fn * z + 3.20092790091004902806E-8;
                    fn = fn * z + 9.41779576128512936592E-11;
                    fn = fn * z + 9.70507110881952024631E-14;

                    fd = 1.00000000000000000000E0;
                    fd = fd * z + 9.17463611873684053703E-1;
                    fd = fd * z + 1.78685545332074536321E-1;
                    fd = fd * z + 1.22253594771971293032E-2;
                    fd = fd * z + 3.58696481881851580297E-4;
                    fd = fd * z + 4.92435064317881464393E-6;
                    fd = fd * z + 3.21956939101046018377E-8;
                    fd = fd * z + 9.43720590350276732376E-11;
                    fd = fd * z + 9.70507110881952025725E-14;
                    f = fn / (x * fd);

                    gn = 6.97359953443276214934E-1;
                    gn = gn * z + 3.30410979305632063225E-1;
                    gn = gn * z + 3.84878767649974295920E-2;
                    gn = gn * z + 1.71718239052347903558E-3;
                    gn = gn * z + 3.48941165502279436777E-5;
                    gn = gn * z + 3.47131167084116673800E-7;
                    gn = gn * z + 1.70404452782044526189E-9;
                    gn = gn * z + 3.85945925430276600453E-12;
                    gn = gn * z + 3.14040098946363334640E-15;

                    gd = 1.00000000000000000000E0;
                    gd = gd * z + 1.68548898811011640017E0;
                    gd = gd * z + 4.87852258695304967486E-1;
                    gd = gd * z + 4.67913194259625806320E-2;
                    gd = gd * z + 1.90284426674399523638E-3;
                    gd = gd * z + 3.68475504442561108162E-5;
                    gd = gd * z + 3.57043223443740838771E-7;
                    gd = gd * z + 1.72693748966316146736E-9;
                    gd = gd * z + 3.87830166023954706752E-12;
                    gd = gd * z + 3.14040098946363335242E-15;
                    g = z * gn / gd;
                }
                si = pi_2 - f * c - g * s;
                if(sg != 0) si = -si;
                ci = f * s - g * c;
            }


            /*************************************************************************
            Hyperbolic sine and cosine integrals

            Approximates the integrals

                                       x
                                       -
                                      | |   cosh t - 1
              Chi(x) = eul + ln x +   |    -----------  dt,
                                    | |          t
                                     -
                                     0

                          x
                          -
                         | |  sinh t
              Shi(x) =   |    ------  dt
                       | |       t
                        -
                        0

            where eul = 0.57721566490153286061 is Euler's constant.
            The integrals are evaluated by power series for x < 8
            and by Chebyshev expansions for x between 8 and 88.
            For large x, both functions approach exp(x)/2x.
            Arguments greater than 88 in magnitude return MAXNUM.


            ACCURACY:

            Test interval 0 to 88.
                                 Relative error:
            arithmetic   function  # trials      peak         rms
               IEEE         Shi      30000       6.9e-16     1.6e-16
                   Absolute error, except relative when |Chi| > 1:
               IEEE         Chi      30000       8.4e-16     1.4e-16

            Cephes Math Library Release 2.8:  June, 2000
            Copyright 1984, 1987, 2000 by Stephen L. Moshier
            *************************************************************************/

            [System.Diagnostics.DebuggerStepThrough]
            public static void HyperbolicSineCosineIntegrals(double x, out double shi, out double chi)
            {
                int sg;
                if(x >= 0) sg = 0;
                else
                {
                    sg = -1;
                    x = -x;
                }

                if(Math.Abs(x) < Eps)
                {
                    shi = 0;
                    chi = -c_MaxRealNumber;
                    return;
                }

                double k;
                double c;
                double s;
                double a;
                if(x < 8)
                {
                    var z = x * x;
                    a = 1;
                    s = 1;
                    c = 0;
                    k = 2;
                    do
                    {
                        a *= z / k;
                        c += a / k;
                        k++;
                        a /= k;
                        s += a / k;
                        k++;
                    } while(Math.Abs(a / s) >= Eps);
                    s *= x;
                }
                else
                {
                    double b1;
                    double b0;
                    double b2;
                    if(x < 18)
                    {
                        a = (576 / x - 52) / 10;
                        k = Math.Exp(x) / x;
                        b0 = 1.83889230173399459482E-17;
                        b1 = 0;
                        ChebIterationShiChi(a, -9.55485532279655569575E-17, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 2.04326105980879882648E-16, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 1.09896949074905343022E-15, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -1.31313534344092599234E-14, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 5.93976226264314278932E-14, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -3.47197010497749154755E-14, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -1.40059764613117131000E-12, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 9.49044626224223543299E-12, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -1.61596181145435454033E-11, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -1.77899784436430310321E-10, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 1.35455469767246947469E-9, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -1.03257121792819495123E-9, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -3.56699611114982536845E-8, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 1.44818877384267342057E-7, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 7.82018215184051295296E-7, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -5.39919118403805073710E-6, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -3.12458202168959833422E-5, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 8.90136741950727517826E-5, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 2.02558474743846862168E-3, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 2.96064440855633256972E-2, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 1.11847751047257036625E0, ref b0, ref b1, out b2);
                        s = k * .5 * (b0 - b2);

                        b0 = -8.12435385225864036372E-18;
                        b1 = 0;
                        ChebIterationShiChi(a, 2.17586413290339214377E-17, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 5.22624394924072204667E-17, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -9.48812110591690559363E-16, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 5.35546311647465209166E-15, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -1.21009970113732918701E-14, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -6.00865178553447437951E-14, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 7.16339649156028587775E-13, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -2.93496072607599856104E-12, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -1.40359438136491256904E-12, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 8.76302288609054966081E-11, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -4.40092476213282340617E-10, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -1.87992075640569295479E-10, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 1.31458150989474594064E-8, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -4.75513930924765465590E-8, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -2.21775018801848880741E-7, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 1.94635531373272490962E-6, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 4.33505889257316408893E-6, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -6.13387001076494349496E-5, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, -3.13085477492997465138E-4, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 4.97164789823116062801E-4, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 2.64347496031374526641E-2, ref b0, ref b1, out b2);
                        ChebIterationShiChi(a, 1.11446150876699213025E0, ref b0, ref b1, out b2);
                        c = k * .5 * (b0 - b2);
                    }
                    else
                    {
                        if(x <= 88)
                        {
                            a = (6336 / x - 212) / 70;
                            k = Math.Exp(x) / x;
                            b0 = -1.05311574154850938805E-17;
                            b1 = 0;
                            ChebIterationShiChi(a, 2.62446095596355225821E-17, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 8.82090135625368160657E-17, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -3.38459811878103047136E-16, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -8.30608026366935789136E-16, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 3.93397875437050071776E-15, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.01765565969729044505E-14, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -4.21128170307640802703E-14, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -1.60818204519802480035E-13, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 3.34714954175994481761E-13, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 2.72600352129153073807E-12, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.66894954752839083608E-12, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -3.49278141024730899554E-11, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -1.58580661666482709598E-10, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -1.79289437183355633342E-10, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.76281629144264523277E-9, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.69050228879421288846E-8, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.25391771228487041649E-7, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.16229947068677338732E-6, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.61038260117376323993E-5, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 3.49810375601053973070E-4, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.28478065259647610779E-2, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.03665722588798326712E0, ref b0, ref b1, out b2);
                            s = k * .5 * (b0 - b2);

                            b0 = 8.06913408255155572081E-18;
                            b1 = 0;
                            ChebIterationShiChi(a, -2.08074168180148170312E-17, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -5.98111329658272336816E-17, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 2.68533951085945765591E-16, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 4.52313941698904694774E-16, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -3.10734917335299464535E-15, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -4.42823207332531972288E-15, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 3.49639695410806959872E-14, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 6.63406731718911586609E-14, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -3.71902448093119218395E-13, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -1.27135418132338309016E-12, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 2.74851141935315395333E-12, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 2.33781843985453438400E-11, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 2.71436006377612442764E-11, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -2.56600180000355990529E-10, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -1.61021375163803438552E-9, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -4.72543064876271773512E-9, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, -3.00095178028681682282E-9, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 7.79387474390914922337E-8, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.06942765566401507066E-6, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.59503164802313196374E-5, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 3.49592575153777996871E-4, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.28475387530065247392E-2, ref b0, ref b1, out b2);
                            ChebIterationShiChi(a, 1.03665693917934275131E0, ref b0, ref b1, out b2);
                            c = k * .5 * (b0 - b2);
                        }
                        else
                        {
                            shi = sg == 0 ? c_MaxRealNumber : -c_MaxRealNumber;
                            chi = c_MaxRealNumber;
                            return;
                        }
                    }
                }
                shi = sg == 0 ? s : -s;
                chi = 0.57721566490153286061 + Math.Log(x) + c;
            }


            private static void ChebIterationShiChi(double x, double c, ref double b0, ref double b1, out double b2)
            {
                b2 = b1;
                b1 = b0;
                b0 = x * b1 - b2 + c;
            }
        }

        /// <summary>(1+sqrt(5))/2</summary>
        private const double c_FibonachiPhi = 1.618033988749895;
        /// <summary>1/sqrt(5)</summary>
        private const double c_InvSqrt5 = 0.4472135954999579;

        [System.Diagnostics.DebuggerStepThrough]
        public static int Fibonachi(int n)
        {
            var phi_n = c_FibonachiPhi.Power(n);
            return (int)(c_InvSqrt5 * (phi_n - Math.Cos(pi * n) / phi_n));
        }

        [System.Diagnostics.DebuggerStepThrough]
        public static int Fibonachi(Complex z)
        {
            var phi_z = c_FibonachiPhi.Power(z);
            return (int)(c_InvSqrt5 * (phi_z - Complex.Cos(pi * z) / phi_z));
        }

        /// <summary>Биномиальный коэффициент</summary>
        /// <param name="n">Показатель степени бинома</param>
        /// <param name="k">Индекс коэффициента</param>
        /// <returns>Биномиальный коэффициент (n, k)</returns>
        public static int BinomialCoefficient(int n, int k)
        {
            if(n >= 0 && k > n || k < 0) return 0;
            if(n < 0) return BinomialCoefficient(-n + k - 1, k) * (k.IsOdd() ? -1 : 1);
            return (int)(n.Factorial() / (k.Factorial() * (n - k).Factorial()));
        }

        /// <summary>Символ Кронекера δ(i,j) = 1 - если i = j, и = 0 - если i ≠ j</summary>
        /// <param name="i">Индекс</param><param name="j">Индекс</param>
        /// <returns>1 - если i = j, и = 0 - если i ≠ j</returns>
        [System.Diagnostics.DebuggerStepThrough, Pure, Axiom]
        public static int KroneckerDelta(int i, int j)
        {
            Contract.Ensures(Contract.Result<int>() == (i == j ? 1 : 0));
            return i == j ? 1 : 0;
        }

    }
}